Test Multivariati e Controllo della Molteplicità tramite Test di Permutazione
livio finos
12/11/2014
Package flip
su CRAN ma soprattutto su github
To install this github version type (in R):
library(devtools)
install_github("flip", "livioivil")
(if devtools is not installed yet: install.packages("devtools") )
seeds data
Leggo i dati seeds
library(flip)
data(seeds)
seeds## grp x y
## 1 0 NA NA
## 2 0 NA NA
## 3 0 NA NA
## 4 0 NA NA
## 5 0 NA NA
## 6 0 NA NA
## 7 0 NA NA
## 8 0 NA NA
## 9 0 6.03 12.54
## 10 0 4.20 14.81
## 11 0 4.49 16.71
## 12 0 2.00 7.53
## 13 0 2.84 7.02
## 14 0 3.88 8.09
## 15 0 2.04 5.76
## 16 0 5.48 18.01
## 17 0 2.31 8.81
## 18 0 1.90 8.17
## 19 0 1.75 6.62
## 20 0 3.02 7.69
## 21 1 NA NA
## 22 1 NA NA
## 23 1 NA NA
## 24 1 3.31 18.49
## 25 1 6.56 19.20
## 26 1 3.16 9.85
## 27 1 4.07 15.83
## 28 1 2.09 6.16
## 29 1 6.72 17.58
## 30 1 3.93 19.29
## 31 1 2.56 10.77
## 32 1 8.30 18.31
## 33 1 4.21 10.56
## 34 1 1.86 9.48
## 35 1 3.09 12.54
## 36 1 5.09 18.35
## 37 1 4.08 11.84
## 38 1 3.63 11.44
## 39 1 2.61 7.66
## 40 1 5.21 12.00# leave out the NAs
dati=seeds[!is.na(seeds$x),]Analisi descrittive:
#funzione che calcola qualche statistica descrittiva
stats=function(v){
sds=apply(v,2,sd)
sds=round(sds,3)
sds=paste("sd :",sds,sep="")
rbind(summary(v),sds)}
#la applico separatamente nei due campioni
by(dati[,2:3],dati$grp,stats)## dati$grp: 0
## x y
## "Min. :1.750 " "Min. : 5.760 "
## "1st Qu.:2.030 " "1st Qu.: 7.402 "
## "Median :2.930 " "Median : 8.130 "
## "Mean :3.328 " "Mean :10.147 "
## "3rd Qu.:4.272 " "3rd Qu.:13.107 "
## "Max. :6.030 " "Max. :18.010 "
## sds "sd :1.467" "sd :4.228"
## --------------------------------------------------------
## dati$grp: 1
## x y
## "Min. :1.860 " "Min. : 6.16 "
## "1st Qu.:3.090 " "1st Qu.:10.56 "
## "Median :3.930 " "Median :12.00 "
## "Mean :4.146 " "Mean :13.49 "
## "3rd Qu.:5.090 " "3rd Qu.:18.31 "
## "Max. :8.300 " "Max. :19.29 "
## sds "sd :1.753" "sd :4.354"box-plot
par(mfrow=c(1,2))
boxplot(seeds$x~seeds$grp,main="x",col=c(0,2))
boxplot(seeds$y~seeds$grp,main="y",col=c(0,2))
Verifica di ipotesi
Confronto:
- H0(z): F(z|grp=0)=F(z|grp=1)
- H1(z): F(z|grp=0)!=F(z|grp=1)
per entrambe le variabili: z=x,y
res=flip(x+y~grp,data=dati,perms=5000)## summary(res)## Call:
## flip(Y = x + y ~ grp, data = dati, perms = 5000)
## 4999 permutations.
## Test Stat tail p-value sig.
## x t 1.320 >< 0.1920
## y t 2.061 >< 0.0532In questo caso è ragionevole fare ipotesi alternative unidirezionali
H1(z): F(z|grp=0)<F(z|grp=1)
res=flip(x+y~grp,data=dati,perms=5000,tail=1)## summary(res)## Call:
## flip(Y = x + y ~ grp, data = dati, tail = 1, perms = 5000)
## 4999 permutations.
## Test Stat tail p-value sig.
## x t 1.320 > 0.0950
## y t 2.061 > 0.0202 *Visualizziamo gli istogrammi delle distribuzioni delle statistiche test di permutazione (calcolate sull’orbita dei dati osservati).
##set the following if you get an error (mostly using Rstudio)
#par(mar=c(1,1,1,1))
hist(res)
E’ però fondamentale studiare la distribuzione congiunta delle due statistiche test:
plot(dati$x,dati$y,pch=21,bg=dati$grp+1)
#medie di gruppo:
XYbars=rbind(colMeans(dati[dati$grp==0,2:3]),colMeans(dati[dati$grp==1,2:3]))
points(XYbars[,1],XYbars[,2],pch=21,bg=c(1:2),cex=2)
abline(h=mean(dati$y))
abline(v=mean(dati$x))
La dipendenza tra le due variabili induce una dipendenza anche nella distribuzione congiunta delle statistiche test
plot(res)
abline(h=0)
abline(v=0)
E di conseguenza anche nei corrispettivi valori lambda (i p-value)
#calcolo i lambda per ogni permutazione, per ogni variabile
res@permP=t2p(res)
#scatter plot:
plot(res@permP[,1],res@permP[,2],col="#F98400",bg="#F2AD00",pch=21)
points(res@permP[1,1],res@permP[1,2],bg="#00A08A",col="#F2AD00",cex=2,pch=21)
abline(v=p.value(res)[1])
abline(h=p.value(res)[2])
ESERCIZIO
è opportuno ricordare che la dipendenza tra statistiche test nella distribuzione congiunta non altera le distribuzioni marginali (cioè univariate) delle statistiche test.
Farte una verifica empirica con comandi del tipo hist(res@permP[,"x"]) o meglio plot(ecdf(res@permP[,"x"]))
Combinazione Non Parametrica di test
Confronto le ipotesi:
- H0: H0(x) AND H0(y)
- H1: H1(x) OR H1(y)
Funzione di Fisher
(res.fish=npc(res,comb.funct = "F")) #sono implementate molte altre funzioni di combinazione e opzioni, si veda ?npc## comb.funct nVar Stat p-value
## V1 Fisher 2 6.256 0.0374Visualizziamo
# calcolo il vettore dei p-values sulla base del vettore di statistiche test
res.fish@permP=t2p(res.fish)
#vettore logico che individua tutti i p-value (ottenuti da permutazioni casuali) che sono <= al p-value combinato osservato
lower=res.fish@permP<=res.fish@permP[1]
#plot
plot(res@permP[,1],res@permP[,2],col="#F2AD00",bg=c("#F98400","#FF0000")[lower+1],pch=21,main="p-values",asp=1)
points(res@permP[1,1],res@permP[1,2],bg="#00A08A",col="#F2AD00",cex=2,pch=21)
#la statistica test osservata è pari a -6.19138 = -log(p-value x) -log(p-value y)
# ( salvato in res.fish@permP[1] )
#la seguente iperbole è il luogo dei punti con la stessa statistica test:
curve(exp(-6.19138)/x,ylim=c(0,1),from=0.000001,to=1,add=TRUE)
Visualizziamo i corrispettivi punti nello spazio delle statistiche test
plot(res,bg=c("#F98400","#FF0000")[(lower)+1],asp=1)
Funzione di Tippett (min-p)
(res.tip=npc(res,comb.funct = "minP")) ## comb.funct nVar Stat p-value
## V1 minP 2 49.5 0.0342Visualizziamo
res.tip@permP=t2p(res.tip)
lower=res.tip@permP<=res.tip@permP[1]
plot(res@permP[,1],res@permP[,2],col="#F2AD00",bg=c("#F98400","#FF0000")[lower+1],pch=21,main="p-values",asp=1)
points(res@permP[1,1],res@permP[1,2],bg="#00A08A",col="#F2AD00",cex=2,pch=21)
abline(h=min(p.value(res)))
abline(v=min(p.value(res)))
Controllo della molteplicità
Errore di I tipo per variabile x : P(p(x)<=alpha|H0)
Nei test esatti questa probabilità è <=alpha
Altrettanto vale per y.
FamilyWise Error Rate (FWER) - Errore di I tipo multiplo
Qual è a probabilità di commettere ALMENO un errore facendo i due test?
P(p(x)<=alpha OR p(y)<=alpha|H0)=?
NON conosciamo la vera configurazione di H0
Non sappiamo cioè se:
- H0(x) AND H0(y) : si possono fare 0, 1 o 2 errori
- H0(x) AND H1(y) : si possono fare 0 o 1 errore
- H1(x) AND H0(y) : si possono fare 0 o 1 errore
- H1(x) AND H1(y) : nessun errore possibile
e quindi 1. H0(x) AND H0(y) : test combinato 2. H0(x) AND H1(y) : test su x (univariato) 3. H1(x) AND H0(y) : test su y (univariato) 4. H1(x) AND H1(y) : nessun test
Bonferroni
- Per controllare FWER a livello alpha, controllo i test univariati a livello alpha/2 (2= numero di test).
- Infatti: P(p(x)<=alpha/2 OR p(y)<=alpha/2|H0)<= P(p(x)<=alpha/2) + P(p(y)<=alpha/2|H0) <=alpha/2 +alpha/2 =alpha
Combinazione di Fisher
plot(res@permP[,1],res@permP[,2],col="#F2AD00",bg="#F98400",pch=21,asp=1)
points(res@permP[1,1],res@permP[1,2],bg="#00A08A",col="#F2AD00",cex=2,pch=21)
# p-value combinato : (caso I)
curve(exp(-6.19138)/x,ylim=c(0,1),from=0.00001,to=1,add=TRUE,col="red",lwd=2)
# p-value per H0(x): (caso II)
abline(v=p.value(res)[1],col="green",lwd=2)
# p-value per H0(y): (caso III)
abline(h=p.value(res)[2],col="blue",lwd=2)
legend("topright",legend=c("H0(x,y) (Fisher)","H0(x)","H0(y)"),col=2:4,bty="n",lwd=2)
Comb. di Tippett (min-p)
plot(res@permP[,1],res@permP[,2],col="#F2AD00",bg="#F98400",pch=21,asp=1)
points(res@permP[1,1],res@permP[1,2],bg="#00A08A",col="#F2AD00",cex=2,pch=21)
# p-value combinato : (caso I)
abline(v=min(p.value(res)),col="red",lwd=2)
abline(h=min(p.value(res)),col="red",lwd=2)
# p-value per H0(x): (caso II)
abline(v=p.value(res)[1],col="green",lwd=2)
# p-value per H0(y): (caso III)
abline(h=p.value(res)[2],col="blue",lwd=2,lty=2)
legend("topright",legend=c("H0(x,y) (Tippett)","H0(x)","H0(y)"),col=2:4,bty="n",lwd=2)
Correzione dei p-value via Fisher (lento quando ci sono molte variabili)
(res=flip.adjust(res,"F"))## Test Stat tail p-value Adjust:Fisher
## x t 1.320 > 0.0950 0.0950
## y t 2.061 > 0.0202 0.0374Correzione dei p-value via minP
(res=flip.adjust(res,"minP"))## Test Stat tail p-value Adjust:Fisher Adjust:minP
## x t 1.320 > 0.0950 0.0950 0.0950
## y t 2.061 > 0.0202 0.0374 0.0342Regioni di rifiuto
Combinazione di Fisher
Combinazione di Tippett (min-p)
